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, Explicit expressions for the transport coefficients of a recently introduced stochastic model for 

' simulating fluctuating fluid dynamics are derived in three dimensions by means of Green-Kubo 
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, simulation data. Two collision rules are considered and their computational efficiency is compared. 
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I. INTRODUCTION 



I In a series of recent papers a discrete-time projection operator technique was used to derive Green-Kubo 

^ ■ relations for the transport coefficients of a new stochastic model — which we will call Stochastic Rotation Dynamics — 
^ ; ' for simulating fluctuating fluid dynamics 4,'F|. Explicit expressions for transport coefficients in two dimensions were 
d derived, and it was shown how random shifts of the collision environment could be used to ensure Galilean invariance 
^ for arbitrary Mach number and temperature. In this paper, we extend our analytical and numerical analysis to 
I three dimensions and consider two distinct collision rules. Expressions for the transport coefficients are derived and 
' ^ compared with simulation results. No assumptions are made regarding molecular chaos, and the correlations which 
can develop at small mean free path are explicitly accounted for. The only approximation we make is to neglect 
fluctuations in the number of particles in the cells which are used to define the collision environment. This amounts 
to neglecting terms of order e~*^, where M is the average number of particles in a cell, and is therefore justified in 
all practical calculations, where M > 5. 

In the Stochastic Rotation Dynamics (SRD) algorithm, the fluid is modeled by particles whose position coordinates, 
ri{t), and velocities, Vi(t), are continuous variables. The system is coarse-grained into cells of a regular lattice, and 
Qs^ there is no restriction on the number of particles in a cell. The evolution of the system consists of two steps: streaming 
, and collision. In the streaming step, all particles are simultaneously propagated a distance v^r, where r is the value 
' of the discretized time step. For the collision step, particles are sorted into cells, and they interact only with members 
' of their own cell. Typically, the simplest cell construction consisting of a hypercubic grid with mesh size a is used. 
As discussed in Refs. 0] and 0, a random shift of the particle coordinates before the collision step is required to 
ensure Galilean invariance. In our implementation of this procedure all particles are shifted by the same random 
vector with components in the interval [—a/2, a/2] before the collision step. Particles are then shifted back to their 
original positions after the collision. If we denote the cell coordinate of the shifted particle i by , the dynamics is 
I , summarized by 

■ r,{t + T) = r,(t)+rv,(i) (1) 

Q : v,(i + r) = uie^t + r)] + u:[e,{t + r)] • {v,:(i) - u[|^(i + r)]}, (2) 



where denotes a stochastic rotation matrix, and u{^^) = -p- J2kei= ^® mean velocity of the particles in 

cell All particles in the cell are subject to the same rotation, but the rotations in different cells are statistically 
independent. There is a great deal of freedom in how the rotation step is implemented, and any stochastic rotation 



^ matrix consistent with detailed balance can be used. The dynamics of the SRD algorithm is explicitly constructed 
- - ' to conserve mass, momentum, and energy, and the collision process is the simplest consistent with these conservation 
laws. The algorithm is Galilean invariant, there is an _ff-theorem, and it yields the correct hydrodynamics equations 
with an ideal gas equation of state [^Q. SRD has been used to study flow around solid objects in both two and 
thr ee p dimensions, dilute polymer solutions B, binary mixtures 10, 11], amphiphilic mixtures 12, 13, 14], coUoids 
[THIIGI] . and cluster structure and dynamics [l^. 

In two dimensions, the stochastic rotation matrix, a;, is typically taken to be a rotation by an angle zta, with 
probability 1/2. Analytic expressions for the transport coefflcients in this case were derived in Refs. and shown 

to be in excellent agreement with simulation results. In three dimensions, one collision rule that has been discussed 
in the literature 0,1a, consists of rotations by an angle a about a randomly chosen direction. All orientations of 
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the random axis occur with equal probabihty. Note that rotations by an angle —a do not need to be considered, since 
this amounts to a rotation by an angle a about an axis with the opposite orientation. The viscosity of this model 
has been measured using a Poiseuille flow geometry in Ref. ■ Analytical expressions for the transport coefficients 
in this case are only available in the limit of large mean free path, A/a — oo, and for one rotation angle, a = 90° 5]. 
In the following, we will refer to this collision rule as Model A. Another, computationally simpler collision rule, which 
we will refer to as Model B, involves rotations about one of three orthogonal rotation axes. In the implementation 
considered here, we take these to be x-, y- and z-axes of a cartesian coordinate system. At each collision step one of 
these three axes is chosen at random, and a rotation by an angle ±q; is then performed, where the sign is chosen at 
random. This procedure is fast and easy to implement; furthermore, only six different rotation matrices are needed, 
which are sparse and contain fixed elements of 1, ± sin(Q;), and cos(a). Our simulations have shown that both collision 
rules lead to a rapid relaxation to thermal equilibrium characterized by the Maxwell-Boltzmann velocity distribution. 

The outline of the paper is as follows. In Sec. ^we briefly summarize the Green-Kubo relations for the transport 
coefficients. Sees. Illll and HVl contain detailed descriptions of the two collision rules, as well as analytical and numerical 
calculations of the shear viscosity, thermal diffusivity and self-diffusion coefficient at both large and small mean free 
path. The work is summarized in Sec. 



II. HYDRODYNAMICS 



The transport coefficients of a simple liquid are the kinematic shear and bulk viscosities, v and 7, and the thermal 
diffusivity coefficient, Dt- Explicit expressions for the asymptotic (long-time limit) shear and bulk viscosities and 
thermal diffusivity of the SRD algorithm were derived in Ref. [2| using a projection operator technique. In particular, 
it was shown that the kinematic viscosities can be expressed in terms of the reduced fluxes in fc-space as 

. + + .i^i = ^ f; V.«(t.o)iA«(t,o). (3) 

while the thermal diffusivity is given by 

Dt^^-j^^Y. {Id+2{kO)\h+2{kt)) , (4) 

where d is the spatial dimension, Cp = {d + 2)^^/2 is the specific heat per particle at constant pressure, and the 
prime on the sum indicates that the t = term has the relative weight 1/2. Here and in the following we have set 
the particle mass equal to one. The thermal conductivity, k, is related to by 

K = pCpDr- (5) 

The reduced fluxes in Eqs. ||2J) and Q are (see Refs. 0,13 for details) 

/i+^k, = ^ E ( - Mm ■ Ae.(o + Av,p{t)k ■ A^m + ^«'(*)) > (6) 

for p = 1 , . . . , d, and 

Id+2{k,t) = ^E('-[(t;f(<)/2-c„r)k.A|,(i) 

i ^ 

+ ^AvUm-Atm+rkBTk-^t)^, (7) 

where — dks/^ is the specific heat per particle at constant volume of an ideal gas, Av'j = Vj{t + t) — v'j{t), 
A|j(i) = |.(t + r) - ^i(i), and A|^(i) ^ ^.{t + t) - |- (i + r), where |j(t + t) is the cell cooridinate of particle i at 
time t + T and £,l{t + r) is the corresponding shifted particle cell coordinate. Since Ari(t) — T\-j{t), Ii{k,t) = to 
this order in k. 
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1. Shear viscosity 

In three dimensions, the shear viscosity can be obtained if, for example, one takes k in the z-direction and [3 — e — 1, 
in the Green-Kubo relation, Eq. ij^l, so that 

oc ' 

'^ = ]^E {hUMHi.t)). (8) 

There are two contributions to the reduced fluxes, namely kinetic and rotational, so that 

/2(z,t)=/|"(z,i)+/2™*(z,t) (9) 

where 

/|™(z,t) = -l^i;„WAe..W (10) 

and 

/r'(z,i) = --5]A^;„(0A^f,(t). (11) 

i 

Contributions to /|™ come from the streaming step, whereas the collisions and shifts contribute to Z™*- There are 
corresponding kinetic, rotational, and mixed contributions to the shear viscosity. 

2. Thermal diffusivity 
Similarly, setting d = 3 and taking k in the z-direction in Eq. one has 

oo ' 

i^T^^r^g (4(z,0)|/.(z,t)) . (12) 

Again, the reduced flux can be divided into the kinetic and rotational contributions, so that 

/5(z,i) =/f-(z,<)+/5™*(z,i), (13) 



where 



and 



""^^^^ ' /^i,S)+rkBTv,,{t)\ (14) 



ir\i,t) = -W\^v]m^tM- (15) 



III. MODEL A: ROTATION AROUND A RANDOM AXIS 

As discussed in the introduction, one choice of collision rule is a rotation by an angle a about a randomly chosen 
axis (see Fig. This collision rule has been used in a recent study of Poiseuille flow and flow around a spherical 
obstacle, and was shown to yield excellent results . Denote the random vector by R; the post-collision velocity of 
a particle at time step t + t can then be written as 

w{t + t) = Uj. {t) + {t) cos(a) + (v'l {t) x R) sin(a) + V|| {t) (16) 

where _L and || denote the components of a vector that are perpendicular and parallel to the random axis R; the 
relative velocity v''(i) = v(t) — U5s(t). 
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The relaxation to thermal equilibrium is characterized by the decay rate of the _ff-function However, a simpler 
procedure is to monitor the relaxation of the fourth moment, S4 = {vf^ + vfy + vf^), of the velocity distribution. 
This was done in Ref. |2| in two dimensions, where it was shown that 5*4 relaxes exponentially to the equilibrium 
value given by the Maxwell-Boltzmann distribution with a relaxation time, r/f, which is essentially temperature 
independent. Furthermore, it was found that is proportional to the average number of particles in a cell, M, and 
depends strongly on the value of the rotation angle a. It diverges approximately as tr ^ a^^ for a — + 0, since there 
are no collisions in this limit, and thermal equilibrium cannot be achieved. As can be seen in Fig. [3 similar behavior 
is observed in three dimensions for both Models A and B. 



A. Large mean free path approximation 

1. Shear viscosity 

For large mean free path. A/a 00, the rotational contributions to the reduced flux, /2°*(z,f), in Eq. © can be 
neglected, so that the shear viscosity can be expressed as 

00 ' 

I' = — — V C„ , (17) 

where 

C„ = (/|*"(z,0)|/|'"(£,nr)) = ^ E {v^AO)^i^z{OM^r)A^,,{nT)). (18) 

As discussed in Ref. except for the t — contribution, Co, it is a good approximation to replace A^^^ by rviz 
when evaluating C„. In the following, we therefore first evaluate v using this approximation, and then discuss the 
required correction term. 

The relevant components of Eq. H16|l can be written as 

Vtxit + r) = u^xit) + c[vi.x{t) ~ u^x{t) -^{vif3{t) - u^,3(t))Ri3Rx] (19) 

+ s[{Viy{t) - U^y{t))Rz - {Viz{t) - U^z{t))Ry\ 

+ X! ("^/^ {t) - u^fj{t))RfjRx, 
Viz (t + T) = U^z (t) + c[Viz it) - u^zit) ~ Y.^v,p{t) - u^p{t))RpRz] (20) 

+ s[{Vix{t) - U^x{t))Ry - {Viy{t) - U^y{t))Rx] 

+ ^{vip{t) - u^i3{t))RpRz, 

13 

where c = cos(a), s = sin(Q;), — '^kee, '^''^ ^^'^ sum runs over all particles in the cell occupied by particle i 
at t = nr. Defining 

A„ ^ {vix{0)viz{0)vjx{nT)vjz{nT)) , (21) 

we have 

^0 = ^{vixVizVjxVjz) = ^{vixVtzVixViz) = N(kBTY, (22) 
ij i 

so that there are only contributions from j — i. The second term in the series is 

Al ='^{vixVizVjxiT)VjziT)), (23) 
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where Vjxij) and Vjz{T) are given by Eqs. H19|l and (|20|l . respectively. There are both diagonal (j = i) and ofF-diagonal 
(j 7^ *) contributions to Ai. Making use of the following averages over the random vector R, 



and 



the diagonal contribution is found to be 



{RD = 1/3 



(i?|i?|,) = 1/15 + 2/15 



{VixVizVix{T)Viz{T)) = {ksTfCA, 



where 



2c(l 



1 ^ 



M' 



(24) 
(25) 
(26) 

(27) 



The ofF-diagonal contribution comes from particles j which are in the same cell as particle i sXt = Q. This contribution 
is equal to 



where 



VA 



,(6c-l)(c-l) 



15Af2 

Since there are M — 1 off-diagonal contributions, it follows that 



(28) 



(29) 



(30) 



The behavior over longer time intervals can be analyzed in a similar fashion. Consider A2 . Following the arguments 
of the last paragraph, there is a diagonal contribution proportional to C^, and an off-diagonal contribution proportional 
to 2(M— l)?7yiCAi since at each time step, M —1 particles become correlated with particle i, and particle j can become 
correlated with particle i at either of the two time steps. Note, however, there are now additional — higher order — 
contributions which arise, for example, when particle j becomes correlated with particle k which then becomes 
correlated with particle i. It is easy to see that these contributions carry additional factors of 1/M and are thus 
of higher order than the diagonal and direct off-diagonal contributions considered above. However, these higher 
off-diagonal contributions can be summed in the geometric series 



(31) 



so that 




C0s(q:) — C0s(2q;)] 



- 1 



(32) 



As discussed above, there is an additional finite cell size correction to this result. It arises from the fact that the 
substitution A^iy = TViy in Co is not precisely correct. Rather, it can be shown that |2J 



C,^A,+N^-^N{kBTf 
6 T 



1 



for A ^ a. Using this result in (|17|l . the corrected kinematic viscosity is 



knT: 



(l-ig)[2-cos(«)-cos(2a)] 



a 

127 



(33) 



(34) 



6 



" 3 






.10^ 


' 12 


©1 



at a « 104.48°. In contrast, in two dimensions, the minimum is at a = 90° for M — > oo, and 

2d _ '^ksT fa 



Note that although this additional term is in general negligibly small in three dimensions, it can dominate the viscosity 
in two dimensions 3]. In particular, for M — > oo, the viscosity in Model A takes the minimum value 

(35) 



C» - ^ i J • (36) 

In this limit, the finite cell size correction provides the sole contribution to the viscosity in two dimensions. The 
viscosity for Model A is always larger than the viscosity in two dimensions. In order to determine the accuracy of 
(|34|l . we have performed simulations using a system of linear dimension L/a — 32, using r = 1, and M = 5 and 20 
particles per cell. Fig. contains a plot of the normalized correlation function (/2(0)/2(i))/iV (fcsT)^ as a function 
of time for two different collision angles, a — 30° and 150°. As expected, the correlations decay much faster for the 
larger collision angle. The resulting time dependent kinematic viscosity is shown in Fig. and the normalized 
asymptotic value of the viscosity, (^/(rfcsT), is plotted in Fig. 2Ji as a function of a for \/a — 2.309, and M = 5 and 
20, and in Fig. ^Jd for A = 1.02 and M = 20. The agreement between the analytical result and simulation data is 
excellent. Fig. contains a plot of the normalized kinematic viscosity, vt/o?', as a function of (A/a)^ for a = 90° 
and M = 20. Also shown in Fig. \^ are data (•) for the viscosity obtained by fitting the one-dimensional velocity 
profile of forced flow between parallel plates in three dimensions ^ . Again, the agreement between both sets of data 
and theory is excellent. 



2. Thermal diffusivity 

As discussed in the previous section, for large mean free path, the rotational contributions to the thermal diffusivity 
in Eq. (|12|l can be neglected. Furthermore, finite cell size corrections of the type discussed in the last section do not 
occur in the calculation of the thermal diffusivity, so that A^^x can be replaced by rviz ■ The thermal diffusivity can 
therefore be expressed as 



CpNkBT^ 



n=0 



(37) 



where B„ = (/|™(z, 0)|/|'"(z, ?it)) with 



JV 



(38) 



Since 



for any value of n, it can be shown that 



{cpT{cpT~ -^)vizVjz{nT)) = 



(39) 



Bo = ^N{kBTf . 

The next term in the series is 

In Appendix A it is shown using quaternion algebra that 



v'f{T)Viz{T)Viz) . 



(40) 



(41) 



(42) 



where 



3 MI 



2 

1 

M 



73 



15M3 



(43) 



7 



and 



7i 



72 



73 



(1 

If, 



2c), 

-l)(8c^ 



3), 



(44) 
(45) 

(46) 



Using quaternion algebra (see Appendix B), it can be shown in the M —>■ oo hmit that the sum in Eq. H37|l is a 
geometric series. Furthermore, direct calculations in two dimensions and for Model B (see Sec. llVAlt suggests 
that this remains true in general. Assuming this is true here, the diagonal contributions to the thermal diffusivity are 
given by 



B, 



= -NikBTfQ-^', 



(47) 



so that carrying out the sum in Eq, 




75 csc2(a/2) 



2 

ksTr 



2 [-64 + 5 A/ (6 

1 / 2 + cos(a) ' 

2 V 1 - cos(a) 



hAf [-3 + 5Af]) + 8 (8- 
3 /4 
M I 5 



■csc^(a/2) 



'5 [-2 + 



M] M) cos(a)] 



- 1 



(48) 
(49) 



Data for the normalized thermal diffusivity, Dtt/o^T, as a function of (A/a)^ for a = 90° and M — is compared 
with (|48|l in Fig. 03. The agreement is excellent. Fig. contains a plot of the various contributions to the time 
dependent correlation function 2{I^{Q)Ic,{t)) /hN{kBT)'^ , and Fig. shows the corresponding data for the time 
dependent thermal diffusivity for a — 30° (filled symbols) and a — f50° (unfilled symbols). Note that for large 
collision angles, stress correlations decay very rapidly, so that only the first couple of terms in the time series are 
needed. Finally, the normalized thermal diffusivity, Dr/rfesT, is plotted as a function of the collision angle in Fig. 
13 for M — 5 and M = 20. Again, the results are in excellent agreement with theory, ft should be emphasized that 
only the diagonal contributions to Dt have been considered here. Although off-diagonal contributions to the thermal 
diffusivity are generally small, better agreement can be achieved for Af < 10 if they are included. In particular, these 
off-diagonal contributions are 0(1/M^). They have been calculated explicitly in two dimensions in Ref. 0, and for 
Model B in Sec. IIVA2l of this paper. 



3. Self-Diffusion Coefficient 
The self-diffusion constant D of particle i is defined by 

D^\iT^^{[v,{t)-vm?)- (50) 
The position of the particle at time t = nr \s 

n-1 

r,(i) =r,(0)-|-r^v,(fcr), (51) 

so that 

n— 1 n—1 

([r,(t) - r,;(0)]^) = ^ E^^(j^) ' ^^(^^))- (52) 

The sums can be rewritten as 

EE (v.(jr).v,(fcr))=^(v,2(jr))-|-2^ ^ (v,(jr) • v,(fcr)) 
3=0 k=o j=a j=o k=j+i 

n-1 

= ndkeT + 2 ^ j(^,{0) ■ v,(n - j)T)). (53) 
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Expression H53|l can be evaluated using the same approximations as were used to determine the viscosity and thermal 
difFusivity. Setting d — 3 and using Eq. lfTC|) , one gets 



(v,(0)-v,(fcT)) =3fcsT7'= , 

where 

7 = [2 cos(a) + l]/3 - 2[cos(a) - 1]/(3M) 
Substituting Ea. (|^ into Ea. lf^ . one gets 

ri-l 



(54) 
(55) 



D = lim ksTr 



1 1 

2 ^ n 



rj ^ — ^ 



n-j 



1+7 ' 



or, explicitly, as a function of M , 



D 



1 - cos(q;) \ M -I 



M 



(56) 



(57) 



The diffusion coefficient was measured for M — 5 and 20 and A/a = 2.309; the results are shown in Fig. |S1 



B. Shear viscosity at small mean free path approximation 

Simple kinetic arguments can be used to calculate the rotational contribution to the kinematic viscosity . Consider 
a collision cell of linear dimension a and divide the cell by the plane z = h. Since the particle collisions occur in 
a shifted cell coordinate system, they result in a transfer of momentum between neighboring cells in the unshifted 
reference frame. The plane z = h represents a cell boundary in the unshifted frame. Consider now the momentum 
transfer in the z-direction, and assume a homogeneous distribution of particles in the cell. The mean velocities in the 
lower and upper partitions are 



^ Ml 



(58) 



and 



U2 



1 



(59) 



respectively, where Mi = M(a — h)/a and M2 — Mh/a. The x-component of the momentum transfer resulting from 
the collision is 



Ml 



Ap^{h) = ^ [vi^x{t + r) - Vi^x{t)] 



(60) 



Using Eq. (|20|l and averaging over the orientation of the vector R then yields 



Since AIu = ^fiUi + M2U2, 



3 a \ a 



so that and averaging over h — which corresponds to averaging over the random grid shift — one has 

(Apa;) = ^(1 - c)M{u2,x - Ui^x)- 



(61) 



(62) 



(63) 
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Since the dynamic viscosity rj is defined as the ratio of the tangential stress, Pzx, to dux/dz, we have 

^_ (Ap,)/(aV) _ (Ap.)/(aV) 

dux/dz (u2.x - ui,x)/(a/2) ' 

so that the kinematic viscosity, v = rj/ p, is 

^=^[l-cos(a)] (65) 

in the Umit of small mean free path. 

We have measured both the rotational and total contributions to the kinematic viscosity for X/a — 0.0361. The 
results are shown in Fig. |^ As can be seen, multi-particle collisions provide the dominant contribution to the viscosity 
for small mean free path. Furthermore, Eq. H65|l provides an accurate approximation for the viscosity in this regime. 
The systematic deviations for small a are due to kinetic contributions (see Fig. 



IV. MODEL B: ROTATION AROUND ORTHOGONAL AXES 



The second collision rule we will consider involves rotations about one of three orthogonal axes. In the implemen- 
tation considered here, we take these axes to be the x-, y-, and z-axes of a cartesian coordinate system. At each 
collision step, one of these axes is chosen at random, and a rotation by an angle ±a is then performed about this axis. 
The sign of a is chosen with equal probability. Rotations about the a;-, y-, and z-axes are described by the rotation 

matrices 

/100\ / c s\ / c s 0\ 

Mx = \ c s , My= \ 10, Af, = -s c , (66) 
\0 -s c J \-s c J \ 1 J 

where c = cos(a) and s = ±sin(Q;), depending on the sign of a. In the following, we will refer to this collision rule 
as Model B. The rate of approach to thermal equilibrium for this model is almost identical to that of Model A. This 
can be seen in Fig. |3 which shows the angular dependence of TfjjM for two values of A/a, 1.15 (•) and 0.0361 (□). 
As in two dimensions, the relaxation rate is essentially independent of temperature. 

An advantage of Model B is that the analytical calculations are comparatively simple and resemble those for the 
model in two dimensions. However, as will be shown in the following section, there are new finite cell size corrections 
which are unique to this collision rule. As will be shown, they occur because rotations are performed about one of 
the symmetry axes of the cell lattice. 



A. Large mean free path approximation 

1. Shear viscosity 

For large mean free path, we proceed as in Sec. IIII A II In order to determine the shear viscosity in this regime, we 
need to evaluate temporal correlation functions of the type 

An = ^ {vix{0)xiy{0)vjx{nT)vjy{nT)) . (67) 

Aq has the same value as in Model A. For n 0, there are again both diagonal {j = i) and off-diagonal [j ^ 
i) contributions to A„. Using the definition of the rotation matrices, Eq. I|66|) . it is easy to show the diagonal 
contributions to Ai are 

A'^=N{kBTfC„ (68) 



Al^N{kBT)\u 

and 

Al=N{kBT)\Cf-e2), 



(69) 
(70) 
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for rotations around the x, y and z-axes, respectively, where Ci = 1/M + c(l — 1/M) and (2 = s{l — 1/M). Averaging 
over the three rotation axes, it follows that the total diagonal contribution is 



where 



^{A^,+Ay + Ai)^N{kBrfCB, 



Cb = (2Ci + Ci - C2)/3. 



(71) 



(72) 



The off-diagonal contributions, which come from particles j which are in the same cell as particle i aX t — Q, can 
be evaluated in a similar fashion. The result is 



where 7?^ = 2c(c - l)/(M-P), so that 



(73) 



(74) 



The off-diagonal contribution is three times smaller than in two dimensions Note that the leading diagonal 
contribution is 0(1), while that of the off-diagonal contribution is 0(1/M). 

The behavior over longer time intervals can be analyzed in a similar fashion, and as for Model A, one finds 



An 



N{kBTf[CB + {M-l)m 



so that 



kBTr f I + Cb + Vb 
I - Cb -riB 



wb{c, 



(75) 



(76) 



The last term on the right hand side of Eq. H76() is a finite cell size correction. In two dimensions and for Model A, 
WA = a^/(12T) 0. As discussed in Sec. IIII A II it occurs because the substitution Af^j, = rviy in the first term on 
the sum on the right hand side of Eq. (jSJ is not precisely correct. In the present case, however, there are additional 
corrections because the rotation matrices always leave one component of the velocity unaltered. As a result, there are 
contributions to Ci that have projections on Cq. 

Finite cell size correction: In order to simplify the discussion of the finite cell size corrections, the following calculations 
are performed in the limit M ^ 00. In this case, the time evolution equations reduce to 



for rotations around the a;-axis. 



for rotations around the j/-axis, and 



Vix{T 



Vix{T 



Vix{T 



Vix{0) 

CViy{0) + ZSViz{0) 



{Q) + ZSV^y{Q) 
CViy{0) - Z (0) 



(77) 



(78) 



(79) 



for rotations around the z-axis, where as before, c and s are the cosine and the sine of the rotation angle a. z = ±1 
specifies the sign of a. 

When calculating Ci, we have to consider rotations about the three symmetry axes separately. As can be seen 
from Eqs. H79|) . rotations about the z-axis mix both the x- and y-components of the velocity, so that the situation is 
similar to that considered in Sec. II B 4 of Ref. Although the same techniques can be used to evaluate Cf as in 
Ref. 0, we know from the results of that paper that there are no finite cell size corrections in this case. 

The situation is different for rotations about the x- and y-axes. For rotations about the x-axis, one has 



T^C^ = fcsT^(A6,(0)A6,(T)), 



(80) 
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so that (A(,iy{0)A^iy{T)) needs to be evaluated. Using the approach described in Sec. II B 4 of Ref. we have 

('[{n+l)a-yo]/T pbi 
l{na~yo)/T 



r°- °° /•[(n+l)a-i/o]/r rbi 

T^Ci /NksT ^ a dyo / dvy / dvznmw{vx)w{vz) , (81) 

Jo ^ ^—^^ J (na~yn)/T J bn 



n.m— — oo 



where ah velocities are at equal time, so that we have dropped the argument (0). Note that the average over z — ±1 
has already be performed. The limits on the inner integral are 

bo = [(m + n)a - yo - Vy{l + c)t]/{st) (82) 

and 

6i = [(m + n + l)a - ?/o - ^^(1 + c)t]/{st). (83) 

w{vx) is the Boltzmann distribution, 

= exp . (84) 



VinksT i 2kBT 

Eq. H81(l looks very similar to Eqs. (18) and (41) in ^ and can be evaluated in an analogous fashion. We therefore 
only sketch the main steps of the analysis, referring to ^] for details. 
The Poisson sum formula [T^ 

E 9in)^ E / gW^-'""'^ dcj^ (85) 

n— — oo m— — oo ^cxd 

is first used twice to transform the double sum over m and n in Eq. (|81|l . Partial integrations over Vx and are 
then performed. The temperature independent part of the resulting expression can be determined by evaluating the 
rh = n = contribution of the remaining sums. The final result of these calculations is 

2 

T^Cf/NkBT^-^ + 0{kBT). (86) 

For rotations about the y-axis, one has 

T^Cf = kBTcJ2{^C^y(.0)^^^yir))■ (87) 

i 

In this case, Viy{T) = Viy{0), and the calculation of {A£^iy{0)A^iy{T)) can be performed using the methods outlined 
above. The final result is 

T^Cf/NkBT = -c^ + 0{kBT). (88) 
Averaging over all three different rotation axes, it follows that 

T^Ci/NkBT = -^il + c) + OikBT). (89) 
do 

Adding this to the contribution from Co, the final approximation for the finite cell size correction for Model B is 



WBic,a'/T)^ — [l--j. (90) 

Although this result is obtained for M ^ oo, and neglects contributions from C„ with n > 2, it reproduces the 
behavior of the viscosity over rotation angles between 10° and 140° and A/a > 0.5 with an error smaller than 2%. 

For M cx), the off-diagonal contributions to the viscosity vanish, and the kinematic viscosity has a minimum at 
a — 120° for A/a oo. For this value of a, 



1 5/0X2 

6 ^ 72 VaJ 



(91) 
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This is significantly smaller than the minimum value given in Eq. (|35|l for Model A, but still larger than the minimum 
value in two dimensions, Eq. Ij^til) . 

Fig. 1101 contains a plot of the normalized kinematic viscosity, (^/(rfcsT), as a function of a for A/a = 1.15 and 
M = 20. Data for the the kinetic (x) and rotational (□) contributions, as well as the total (•) viscosity, are plotted 
and compared with the theoretical prediction, Eqs. IjTGI) and (|90|l . The agreement is excellent. Note, in particular, 
that the finite cell size contribution to the total viscosity is not negligible, particularly for large rotation angles. In 
Fig. ^] the normalized viscosity v/{TkBT) is plotted as a function of M for a = 90° and A/a = 1.15. Again, the 
agreement between theory and simulation is excellent. Finally, Fig. 1121 shows the normalized total shear viscosity, 
vt/o?, as a function of (A/a)^ for M — 20 and a = 90°. Note in particular that the both the M dependence of 
the viscosity as well as the size of the finite cell size correction — given by the intercept — are accurately described by 
theory. 



2. Thermal diffusivity 

The kinetic part of the reduced flux for the calculation of the thermal diffusivity is given by Eq. (|38|l . where again, 
Bn = (/^^"(z, 0)|/|™(z, tit)). The calculation of the thermal diffusivity simplifies considerably if we utilize relation 
(|39|l and the following relations, 

CpT - <\ ==0 for m = 1, 2, and 3, (92) 



and 



CpT - 



CpT 



vl = 3(fcB^)^ 



Bo is the same as in Model A, namely 



Bo = ^NikeTf. 



Bi (including off-diagonal terms) is 



= -N{kBTf [774 + (Af - 1)75] , 



where 7 is given given in Eq. (|55|l . 

74 = [1 + 2(C? + C2')]/3, 
where C,i and C,2 are defined in the text following Eq. (|70|l . and 

16(l-c)2 



(93) 



(94) 



(95) 



(96) 



(97) 



75 



15 AP 



The coefficients _B„ form a geometrical series, because successive rotations are uncorrelated. This can be seen by 
first performing an average over the rotation angle and then performing the thermal average. In particular. 



Bn = ^N{kBTf [774 + [M - 1)75]" , 



so that the thermal diffusivity is 



Dt = 



ksTr 



1 + 774 + (M - 1)75 
1 - 774 - (M - 1)75 



(99) 



(100) 



Note that the off-diagonal contribution is of order l/A'P; it is therefore less important than for the shear viscosity. 
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3. Self-Diffusion Coefficient 



The diffusion constant can be determined as in Sec. IIII A 21 for Model A. The final result is 

which is the same as for Model A [see Eq. ()57|) ]. It is interesting to note that 

Dt 

-jj 1 forAf ^ oo, (102) 
for both Model A and B as well as in two dimensions. 



B. Small mean free path approximation: shear viscosity 



A detailed calculation of the_shear viscosity in this limit can be performed following the arguments used in Sec. 
IIII A 21 for Model A and in Ref. 3] for two-dimensions. However, for Model B, the following simple argument gives the 
same result. Consider the momentum transfer across a plane perpendicular to the z-axis. Only rotations about the 
X- and y-ax.es produce a nonzero momentum transfer, and since the momentum transfer — and therefore the resulting 
viscosity — from each of these rotations is equal to that calculated two-dimensions 3] , one finds that 

2 

'^3D = i:V2D = -^[1 - cos(q!)]. (103) 
3 18t 

Note that this expression is identical to the one obtained for Model A. Data for the a dependence of the normalized 
viscosity, vt/o^, at A/a — 0.0361 is plotted in Fig. E| Note in particular the importance of kinetic contributions to 
the viscosity for small a, even for this small value of A/a. 



V. SUMMARY 



In this paper we have presented a comprehensive analytical and numerical study of the stochastic rotation dynamics 
model for fiuid dynamics in three dimensions for two collision rules. The first collision rule (Model A) consists of a 
rotation by an angle a about a randomly chosen axis. It was introduced in Refs. 0| and and used in Ref. jH to 
study channel flow and flow about a spherical object. A new, simpler collision rule (Model B), in which collisions 
involve rotations by an angle ±a about one of three orthogonal axes, was also discussed. Calculations involving this 
model are particularly simple, since the rotations about the individual axes are very similar to those in two dimensions. 
In particular, it was possible using this model to calculate the off-diagonal contributions to the thermal diffusivity; 
a similar calculation for Model A was prohibitively tedious. Since both models are comparable with regard to their 
computational efficiency, i.e. relaxation rates, range of viscosities, etc., the simplicity of Model B can have advantages 
in specific applications. 

Discrete time Green-Kubo relations originally derived in Refs. 0] and were used to determine explicit expressions 
for the shear viscosity, the thermal diffusivity, and the self-diffusion constant. The kinetic, collision, and mixed 
contributions to the transport coefficients were analyzed individually, and no assumptions regarding molecular chaos 
were made. This enabled us to determine correlation induced finite cell size corrections to the shear viscosity which 
persist even in the limit of large mean free path. In Ref. |H it was shown that these corrections can, under certain 
circumstances, such as collisions with a = 90° and large particle density, provide the dominant contribution to the 
shear viscosity in two dimensions. In three dimensions, we showed here that corrections of this type, while not entirely 
negligible, are rather small for Model A. However, as discussed in Sec. IIV A II for Model B, where collisions involve 
rotations about one of three previously defined orthogonal axes, there are additional finite cell size corrections which 
make non-negligible contributions to the viscosity for a wide range of densities and rotation angles. It is important 
to note that corrections of this type are only important for the shear viscosity. 

It was also shown how quaternion algebra can be used to simplify calculations of kinetic contributions to the 
transport coefficients. In particular, the appendicies describe the calculation of the thermal diffusivity in Model A 
using quaternions. Finally, simulation results for the viscosity, thermal diffusivity, and the self-diffusion coefficient for 
range of simulation parameters were presented and compared to the analytical approximations. In all cases, agreement 
was excellent; furthermore, the comparisons showed that the finite cell size corrections described above are necessary 
in order to achieve quantitative agreement. 
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APPENDIX A 

The calculation of correlation functions of the reduced fluxes can be simplified by rewriting the time evolution 
equations for the velocities using quaternions. Two arbitrary quaternions, P and Q, are defined by 

P = (p,P) (Al) 
Q = ('Z,Q), (A2) 

where {p,q} are the scalar parts and {P,Q} the corresponding vector parts of the quaternions. If the scalar part 
is zero, the quaternion is an ordinary vector and is called as a "pure" quaternion. The multiplication rule of two 
quaternions is given by 

PQ= (re-P.Q, pQ + gP + P X Q) . (A3) 
It follows that for two pure quaternions, R = (0, R) and $ = (0, S), 

IRm= (0,-|R|2S) . (A4) 

Defining 

¥(i) = (0,v(t)) , (A5) 
U=(0,u,), (A6) 

and 

¥'■ = (0,v'-) = ¥(0) -U, (A7) 
The time evolution equation for the velocities, Eq. (|l(j|l . can be written as 

¥(t) ^ A¥''A* + U . (AS) 

where 

A = (cos(a/2), Rsin(a/2)) . (A9) 



The first term in Eq. (jA8|) corresponds to the rotation of the relative velocity vector around the random axis R. 
Using the multiplication rule given in Eq. ljA3|l . it is easy to see that Eq. (|A8|) is equivalent to Eq. H16|l . Similarly, 
using Eq. (|A4|I it can be shown that 

{VHt)),=-v^{t)v,{t) . (AlO) 
Dropping the index i, Bi given by Eq. (|41() can be written as 

^ = \{v'v'{r)v^ir)v,)-^{v^ir)v.ir)v.), (All) 

" V ' " V ' 

B[ B[' 

or by using Eq. HA10|) . as 

^ = \{{V\r)).v.{vl+vl + vl))-"-^{{VHr)).v.) • (A12) 
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Using the multiplication rule for quaternions, and the fact that AA* = 1, it can be shown that 

¥^(t) = A {vyj?C + U^AV^A + AVAUAVA + UA(V)2A 
+ A(¥'')2A*U + U^ + A¥''AU^+UA¥'^AU . 

Simplifying terms and using energy conservation, 



(A13) 
(A14) 



(A15) 



one obtains 



1 



B[ = -riivl +vl+vl + 2us^{v^{t) - v^) + 2u^y{vy{T) - Vy) + 2us^(w^(t) - vj){vl + + vI)v,{t)v,) (A16) 



{ksTf 



12 



35(1 + 2c) 



2(1 -c) 
M 



20, 

31- 16c + — 2c- 1 



576(1 - cf 
5M3 



and 



{{^l + + 2u^x{vxiT) - Vx) + 2u^y{vy{T) - Vy) + 2u^2{v^iT) ~ Vz))v^{t)v^) (A18) 



12 



^ 10, ^ 16 , 4 , 

^(1 + 2^)+m(1-^)-5M^(^-'^)(1-m) 



which then yields Eq. when substituted into Eq. |IA11|I . 



(A17) 



(A19) 



APPENDIX B 

In the limit Af — > oo, U 0, and Eq. IjASjl can be written as 

¥(t) = A¥A, 

where we have dropped the superscript "r", so that ¥ = ¥(0). The cube of ¥(t) is then simply 

¥^(t) = A¥^A*, 

where 

= (o,-|v|2v) . 



(Bl) 
(B2) 
(B3) 

This means that ¥'^(r) is the rotation of the vector — |vpv around a random axis R. Eqs. (jB2p and (jB3p can be used 
to evaluate the second term in Eq. IjAlip . namely 



El 



CpT 



which can be shown to equal 



2 cos(a) + 1 



(B4) 



(B5) 



Similarly, for t = 2t, 

¥(2r) = A¥(r)A*, (B6) 
where prime denotes a different random vector then in Eq. l)A9p . Using energy conservation and the commutator 

[A,¥] = (0,2sin(a/2) R' X v), (B7) 
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V^{2t) can be written as 



V^(2r) = -|v(2t)|2¥(2t) 

= -\v{t)\^KV{t)K* 

= -|v(r)|2(¥^ + [A',¥])A'* 

= -|v(t)|2¥- |v(r)|2(0,2sin^ R' x v)A'* 



(B8) 

(B9) 
(BIO) 

(Bll) 



so that 



Since 



one gets finally, 



2 cos(a) + 1 



E2 = 



2cos(a) + 1 



£1 = 5(fcsT)2 



2 cos(a) + 1 



(B12) 



(B13) 



(B14) 



so that the terms B'^ form a geometric series. It can also be shown that the B!^ are terms in a geometric series, with 
the same angular dependence. The difference of these two terms is therefore also a geometric series. 
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FIG. 1: Rotation of the vector v'' around a random direction R by the angle a. 
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FIG. 2: The normalized relaxation time, trJM, of the fourth moment of the velocity distribution, 5*4 — + "^ty + I'iz) as 

a function of the rotation angle a for M = 20, where M is the average number of particles per cell. The data for Model A (*) 
were obtained for A/a = 1.15, while the data for Model B correspond to A/a = 1.15 (•) and A/a = 0.0361 (□). Parameters: 
L/a = 32 and r = 1. 
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FIG. 3: (a) Normalized correlation functions {l2{Q)h{i)) /N{kBTf for Model A as a function of time for a = 30° (solid 
symbols) and a = 150° (unfilled symbols). For a = 30°, the kinetic, rotational, and mixed contributions arc indicated by •, ■, 
and <, respectively. For a = 150°, the kinetic, rotational, and mixed contributions are indicated by o, and O, respectively, 
(b) Normalized time dependent kinematic viscosity, v{t) /rkBT . Symbols are the same as in part (a). Parameters: L/a = 32, 
A/a = 2.309, r = 1, and M = 20. The data were obtained by time averaging over 75000 iterations. 
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FIG. 4: Normalized kinematic viscosity, u/rkBT, for Model A as a function of the collision angle a. (a) Data for L/a — 32, 
A/a = 2.309, r = 1, and M = 5 (■) and M = 20 (•). (b) Data for L/a = 32, A/a = 1.02, r = 1, and M = 20. The bullets 
are results obtained using the Green-Kubo relation, and the unfilled boxes (□) are data for the kinematic viscosity obtained 
in Ref. Q by fitting the one-dimensional velocity profile of forced flow between parallel plates. The lines are the theoretical 
prediction, Eq. 1341 . for the corresponding parameter values. The data were obtained by time averaging over 75000 iterations. 
The deviation of the data point □ at a = 30° is due to finite Knudsen number effects. 
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FIG. 5: a) Normalized kinematic viscosity, vt /a^ , and b) thermal diffusivity, Dtt /a^ , for Model A as functions of (A/a)^ for 
collision angle a — 90°. The bullets are data obtained using Green-Kubo relations. The unfilled boxes (□) are data for the 
kinematic viscosity obtained in Ref. |3 by fitting the one-dimensional velocity profile of forced flow between parallel plates. 
The solid line is the theoretical prediction, Eqs. 1341 and 14811 . Parameters: L/a = 32, t = 1, and M = 20. The data were 
obtained by time averaging over 75000 iterations. 
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FIG. 6: (a) Normalized correlation functions 2{h{0)h{t)) /bN(kBTf for Model A as a function of time for a = 30° (filled 
symbols) and a = 150° (unfilled symbols). For a = 30°, the kinetic, rotational, and mixed contributions arc indicated by •, ■, 
and respectively. For a = 150°, the kinetic, rotational, and mixed contributions are indicated by o, and O, respectively, 
(b) Normalized time dependent thermal diffusivity, Drit) /rkBT . Symbols arc the same as in part (a). Parameters: L/a = 32, 
A/a = 2.309, r = 1, and M = 20. The data were obtained by time averaging over 75000 iterations. 
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FIG. 7: Normalized thermal diffusivity, Dt /TknT , for Model A as a function of collision angle a. The lines are the theoretical 
prediction, Eq. (I48II . The data were obtained by time averaging over 75000 iterations. Parameters: L/a = 32, A/a = 2.309, 
T = 1, and M = 5 (o) and M = 20 (•). 
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FIG. 8: Normalized self-diffusion constant, D/rkBT, for Model A as a function of collision angle a. The lines are the theoretical 
prediction, Eq. 1571 . The data were obtained by tme averaging over 75000 iterations. Parameters: L/a = 32, A/a = 2.309, 
r = 1, and M = 5 (■) and M = 20 (•). 
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FIG. 9: Normalized kinematic viscosity for Model A as a function of collision angle a for small mean free path, A/a = 0.0361. 
The plot shows both the rotational (•) and the total (■) contributions to the viscosity. The solid line is the theoretical 
prediction, Eq. 13411 . The data were obtained by time averaging over 75000 iterations. Parameters: L/a — 32 and M — 20. 
The open squares (□) are data for the total kinematic viscosity obtained in Ref. Q. 
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FIG. 10: Various contributions to the normalized shear viscosity, u/rkBT, for Model B as a function of the rotation angle a at 
large mean free path, A/a = 1.15. The symbols are simulation data for the kinetic contribution (x), the rotational contribution 
(□), and the total viscosity (•). The solid line is the theoretical prediction, Eqs. 17611 and 19011 . The viscosity has a minimum 
at a = 120°, as predicted by theory. The data were obtained by time averaging over 40000 iterations. Parameters: L/a = 32, 
r = 1, and M = 20. 
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FIG. 11: Various contributions to the normalized shear viscosity, u/rkBT, for Model B as a function of M for large mean free 
path, A/a = 1.15, and a — 90°. The symbols are simulation data for the kinetic contribution (x), the rotational contribution 
(□), and the total viscosity (•). The solid line is the theoretical prediction, Eqs. 1761 and 1901 1. For this value of A/a, rotational 
contributions to the total viscosity are negligible. 
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FIG. 12: Shear viscosity for Model B as a function of A/a for a = 90° and M = 20. The symbols are simulation data for the 
kinetic contribution (□) and the total viscosity (•). The slope of 0.297 is in excellent agreement with the theoretical prediction, 
0.2895, which follows from Eqs. 1761 and 19011 . Note that for M — » oo, theory predicts a slope of 0.25; 1/M corrections are 
therefore important even for M = 20. Parameters: L/a — 32, r = 1. 
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FIG. 13: Normalized shear viscosity, vt/o^, for Model B as a function of the rotation angle a for small mean free path, 
A/a = 0.0361. The bullets are simulation data and the solid line is the theoretical prediction, Eq. I103L for the rotational 
contribution to the kinematic viscosity. The deviation of the data from the theoretical prediction for a < 30° is due to the 
increasing importance of the kinetic contribution to the viscosity (see Fig. 110^ . The data were obtained by time averaging over 
40000 iterations. Parameters: L/a = 32, t = 1 and M = 20. 



